/* benchmark guru interfaces */
#include <benchmark/fft_bench.h>
#include <Configuration.h>
#include <mpi/communicator.h>
#include <mpi/collectives.h>
#include <Message/OpenMP.h>
#include <Utilities/Timer.h>
#include <bitset>

using namespace APPNAMESPACE;

int main(int argc, char** argv)
{

  mpi::environment env(argc,argv);
  mpi::communicator mycomm;
  OhmmsInfo ohmms("fft_omp",mycomm.rank(),0,1);

  typedef  double real_type;
  typedef  std::complex<double> complex_type;

  ///1D FFT dimension
  int fft_dim=12288;
  ///total number of ffts : 512*16 cores
  int fft_tot=8192;
  ///number of fft blocks per call
  int fft_block=4;
  int niters=10;
  unsigned plan_mode = FFTW_MEASURE;

  std::bitset<MAX_FFT_ENG> test_it(0);
  test_it.set(FFTW_ENG,1);
  vector<string>  lib_names;
  lib_names.push_back("FFTW");
  lib_names.push_back("ESSL");
  lib_names.push_back("MKL");

  //cout.setf(std::ios::scientific, std::ios::floatfield);
  int ic=0;
  while(ic<argc)
  {
    string a(argv[ic]);
    if(a.find("-s")<a.size()) //get fft dimension
      fft_dim=atoi(argv[++ic]);
    else if(a.find("-n")<a.size()) //get how many fft to perform
      fft_tot=atoi(argv[++ic]);
    else if(a.find("-b")<a.size()) //number of simultaneous trans.
      fft_block=atoi(argv[++ic]);
    else if(a.find("-i")<a.size()) 
      niters=atoi(argv[++ic]);
    else if(a.find("disable-fftw")<a.size())
      test_it.set(FFTW_ENG,0);
#if defined(HAVE_LIBESSL)
    else if(a.find("enable-essl")<a.size())
      test_it.set(FFTESSL_ENG,1);
#endif
#if defined(HAVE_LIBMKL)
    else if(a.find("enable-mkl")<a.size())
      test_it.set(FFTMKL_ENG,1);
#endif
    ++ic;
  }

  //fft per node
  fft_tot=fft_tot/mycomm.size();
  //set the block of pencils
  fft_block=std::min(fft_block,fft_tot);

  //allocate the fft engines per thread
  int max_threads=omp_get_max_threads();
  typedef FFT_DFT<real_type,real_type,true,FFTW_ENG> fft_engine_t;
  std::vector<fft_engine_t*> ffts(max_threads);
  //std::vector<ofstream*> fout(max_threads);
#pragma omp parallel 
  {
    vector<real_type> target((fft_dim+2)*fft_block);
    ffts[omp_get_thread_num()]=new fft_engine_t(fft_dim,&target[0],fft_block);
    //char fname[16];
    //sprintf(fname,"test_%i",omp_get_thread_num());
    //fout[omp_get_thread_num()]=new ofstream(fname);
  }

  //loop count for a fft
  int fft_loop=fft_tot/fft_block;


  double mflops_c=(std::log(static_cast<double>(fft_dim))/std::log(2.))*5.*fft_dim*fft_tot*1e-6;
  double mflops_r=mflops_c*0.5;
  app_log() << "#  Testing FFT dim = " << fft_dim 
    << " total ffts = " << fft_tot*mycomm.size()
    << " per node = " << fft_tot
    << " per thread = " << fft_tot/max_threads
    << " fft bundles = " << fft_block
    << " fft loop =" << fft_loop
    << " iterations= " << niters 
    << " MFLOPS(complex)= " << static_cast<int>(mflops_c)
    << endl;

  //allocate the real data
  vector<real_type> in((fft_dim+2)*fft_tot);

  int offset=(fft_dim+2)*fft_block;

  Timer myclock;
  myclock.restart();
  for(int iter=0; iter<niters; ++iter)
  {
#pragma omp parallel for schedule(static,4)
//#pragma omp parallel for 
    for(int i=0; i<fft_loop; ++i)
    {
      ffts[omp_get_thread_num()]->transform_forward(&in[0]+i*offset,fft_block);
      //*fout[omp_get_thread_num()] << i*offset<<endl;
    }
  }

  double wct=myclock.elapsed();
  mpi::reduce(mycomm,wct);

  app_log().setf(ios::floatfield, ios::scientific);
  app_log().precision(3);
  app_log() 
    << " FINAL " << fft_tot*mycomm.size()
    << " " << fft_block
    << " " << fft_loop
    << " " << niters 
    << " WCT " << wct/static_cast<double>(mycomm.size()*niters) << endl;

  for(int i=0; i<ffts.size(); ++i) delete ffts[i];

  mpi::environment::finalize();
  return 0;
}

